function plotTermPremia(outEstim)

mat = 40; %Maturity for Term Premia
setupFilter     = outEstim.setupFilter;
timeIndexPlot   = 1965:10:2015;  
fontSizeValue   = 14;
figureSize      =  [0 0 1300 600+300]; % (0,0) size in x size in y

termPrem        = outEstim.termPrem;
termPremReal.TP = termPrem.RealTP;
termPremReal.MeanTP = termPrem.RealMeanTP;

%% Plot Nominal Term premium 
version = 1;
figure('Name','Nominal Term Premium','NumberTitle','off','Position',figureSize)
hold on
NBER        = outEstim.setupFilter.NBERrec(setupFilter.numInitial+1:end);
yMax        = max([termPrem.TP(mat,setupFilter.numInitial+1:end),setupFilter.acm_tp.tp10y']);
yMin        = min([termPrem.TP(mat,setupFilter.numInitial+1:end),setupFilter.acm_tp.tp10y']);
dates       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
if version == 1
    xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
    nn          = length(dates);
    idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
    reverseaxis = setupFilter.timeIndex(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin*100];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
else
    %For shading above zero
    X  = [dates',fliplr(dates)'];
    Y1 = [100*yMax*NBER' 0*NBER'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);
    
    % For shading below zero
    Y2 = [100*yMin*NBER' 0*NBER'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);
end
hold on
tmp1 = plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),termPrem.TP(mat,setupFilter.numInitial+1:end),'-k','LineWidth',2);
tmp2 = plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
    setupFilter.acm_tp.tp10y(setupFilter.numInitial+1:end),'-k');
legend([tmp1(1), tmp2(1)],{'$M^{M,CS}$','Adrian et al (2013)'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'interpreter','Latex','FontSize',fontSizeValue+2);
hold off

%axis tight
axis([dates(1) dates(end) yMin yMax])
datetick('x','yyyy','keeplimits','keepticks')
xticks(datenum(datetime(timeIndexPlot,1,1)))
xticklabels(num2cell(timeIndexPlot))
set(gca,'FontSize',fontSizeValue)
corrTPnom_ACM = corr(termPrem.TP(mat,setupFilter.numInitial+1:end)',setupFilter.acm_tp.tp10y(setupFilter.numInitial+1:end))
title([num2str(mat/4),'-year term premium: Mean = ',num2str(termPrem.MeanTP(mat,1),3)],'FontSize',fontSizeValue+2)

%% Plot Real Term premium 
version = 1;
figure('Name','Real Term Premium','NumberTitle','off','Position',figureSize)
hold on
yMax        = max(termPremReal.TP(mat,setupFilter.numInitial+1:end));
yMin        = min(termPremReal.TP(mat,setupFilter.numInitial+1:end));
if version == 1
    xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
    nn          = length(dates);
    idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
    reverseaxis = setupFilter.timeIndex(idx);
    index       = 1:nn;
    reversIndex = nn:-1:1;
    grpyat = [xaxis       100*yMax*NBER(index,1);...
        reverseaxis NBER(reversIndex,1)*yMin*100];
    patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
else
    %For shading above zero
    X  = [dates',fliplr(dates)'];
    Y1 = [100*yMax*NBER' 0*NBER'];
    p2 = fill(X,Y1,[0.8 0.8 0.8]);
    set(p2,'edgecolor',[0.8 0.8 0.8]);
    
    % For shading below zero
    Y2 = [100*yMin*NBER' 0*NBER'];
    p3 = fill(X,Y2,[0.8 0.8 0.8]);
    set(p3,'edgecolor',[0.8 0.8 0.8]);
end
plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),termPremReal.TP(mat,setupFilter.numInitial+1:end),'-k')

%axis tight
axis([dates(1) dates(end) yMin yMax])
datetick('x','yyyy','keeplimits','keepticks')
xticks(datenum(datetime(timeIndexPlot,1,1)))
xticklabels(num2cell(timeIndexPlot))
set(gca,'FontSize',fontSizeValue)
title([num2str(mat/4),'-year real term premium: Mean = ',num2str(termPremReal.MeanTP(mat,1),3)],...
    'FontSize',fontSizeValue+2)

% % %% Plot Inflation Risk Term premium 
% % version = 1;
% % figure('Name','Inflation Risk Term Premium','NumberTitle','off','Position',figureSize)
% % hold on
% % yMax        = max(termPrem.lTPxhat(mat,setupFilter.numInitial+1:end));
% % yMin        = min(termPrem.InflTPxhat(mat,setupFilter.numInitial+1:end));
% % if version == 1
% %     xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
% %     nn          = length(dates);
% %     idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
% %     reverseaxis = setupFilter.timeIndex(idx);
% %     index       = 1:nn;
% %     reversIndex = nn:-1:1;
% %     grpyat = [xaxis       100*yMax*NBER(index,1);...
% %         reverseaxis NBER(reversIndex,1)*yMin*100];
% %     patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);    
% % else
% %     %For shading above zero
% %     X  = [dates',fliplr(dates)'];
% %     Y1 = [100*yMax*NBER' 0*NBER'];
% %     p2 = fill(X,Y1,[0.8 0.8 0.8]);
% %     set(p2,'edgecolor',[0.8 0.8 0.8]);
% %     
% %     % For shading below zero
% %     Y2 = [100*yMin*NBER' 0*NBER'];
% %     p3 = fill(X,Y2,[0.8 0.8 0.8]);
% %     set(p3,'edgecolor',[0.8 0.8 0.8]);
% % end
% % plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),termPrem.InflTPxhat(mat,setupFilter.numInitial+1:end),'-k')
% % 
% % %axis tight
% % axis([dates(1) dates(end) yMin yMax])
% % datetick('x','yyyy','keeplimits','keepticks')
% % xticks(datenum(datetime(timeIndexPlot,1,1)))
% % xticklabels(num2cell(timeIndexPlot))
% % set(gca,'FontSize',fontSizeValue)
% % title([num2str(mat/4),'-year inflation risk term premium: Mean = ',num2str(termPrem.meanInflTP(mat,1),3),...
% % ' Std = ',num2str(termPrem.stdInflTP(mat,1),3)],'FontSize',fontSizeValue+2)
% % 
%% Shock Decomposition of Nominal Term Premia
version = 1;
figure('Name','Nominal Term Premium: Decomposition','NumberTitle','off','Position',figureSize)
hold on
NBER        = outEstim.setupFilter.NBERrec(setupFilter.numInitial+1:end);
yMax        = max(max(termPrem.TPdecom(mat,setupFilter.numInitial+1:end,:)));
yMin        = min(min(termPrem.TPdecom(mat,setupFilter.numInitial+1:end,:)));
dates       = setupFilter.timeIndex(setupFilter.numInitial+1:end);

termPremDeMean_TPdecom = termPrem.TPdecom - repmat(mean(termPrem.TPdecom,2),1,size(termPrem.TPdecom,2),1);
totalAbs      = sum(abs(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:)),3);
fracShockAbs  = squeeze(sum(abs(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:)))./sum(totalAbs));
totalRMSE     = sum((termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:).^2),3);
fracShockRMSE = squeeze(sum(termPremDeMean_TPdecom(mat,setupFilter.numInitial+1:end,:).^2)./sum(totalRMSE));

for i=1:outEstim.model.ne
    subplot(2,3,i)
    if version == 1
        xaxis       = setupFilter.timeIndex(setupFilter.numInitial+1:end);
        nn          = length(dates);
        idx         = nn+setupFilter.numInitial:-1:setupFilter.numInitial+1;
        reverseaxis = setupFilter.timeIndex(idx);
        index       = 1:nn;
        reversIndex = nn:-1:1;
        grpyat = [xaxis       100*yMax*NBER(index,1);...
            reverseaxis NBER(reversIndex,1)*yMin*100];
        patch(grpyat(:,1),grpyat(:,2),[0.8 0.8 0.8],'edgecolor',[0.6 0.6 0.6]);
    else
        %For shading above zero
        X  = [dates',fliplr(dates)'];
        Y1 = [100*yMax*NBER' 0*NBER'];
        p2 = fill(X,Y1,[0.8 0.8 0.8]);
        set(p2,'edgecolor',[0.8 0.8 0.8]);
        
        % For shading below zero
        Y2 = [100*yMin*NBER' 0*NBER'];
        p3 = fill(X,Y2,[0.8 0.8 0.8]);
        set(p3,'edgecolor',[0.8 0.8 0.8]);
    end

    hold on
    %plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
    %    termPrem.TP(mat,setupFilter.numInitial+1:end)-termPrem.TPdecom(mat,setupFilter.numInitial+1:end,i),'-k');
    plot(setupFilter.timeIndex(setupFilter.numInitial+1:end),...
        termPrem.TPdecom(mat,setupFilter.numInitial+1:end,i),'-k');
    hold off
    axis([dates(1) dates(end) yMin yMax])
    datetick('x','yyyy','keeplimits','keepticks')
    xticks(datenum(datetime(timeIndexPlot,1,1)))
    xticklabels(num2cell(timeIndexPlot))
    set(gca,'FontSize',fontSizeValue-2)
    title([outEstim.model.labelx{outEstim.model.nx1+i},...
        ': pct($\Delta TP_t$) = ',num2str(round(100*fracShockAbs(i),1))],...
        'Interpreter','Latex','FontSize',fontSizeValue)
    %title([outEstim.model.labelne{i},': pct($\Delta TP_t$) = ',num2str(round(100*fracShockRMSE(i),1))],'Interpreter','Latex','FontSize',fontSizeValue)
end


end